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Abstract 

The paper presents a time-domain method for com- 
putation of sound radiation from aircraft engine 
sources to the far-held. The effects of nonuniform 
how around the aircraft and scattering of sound by 
fuselage and wings are accounted for in the formu- 
lation. Our approach is based on the discretization 
of the inviscid how equations through a collocation 
form of the Discontinuous Galerkin spectral element 
method. An isoparametric representation of the un- 
derlying geometry is used in order to take full advan- 
tage of the spectral accuracy of the method. Large- 
scale computations are made possible by a parallel 
implementation based on message passing. Results 
obtained for radiation from an axisymmetric nacelle 
alone are compared with those obtained when the 
same nacelle is installed in a generic conhguration, 
with and without a wing. 

I. Introduction 

One of the main sources of aircraft noise is the 
engine. Most of the measurements and model- 
ing regarding engine noise have been done for 
an isolated engine [3, 9]. The far field engine 
noise of an aircraft is, however, influenced both 
by the flow around the wings and fuselage and 
by scattering of sound from aircraft surfaces. 
Ideally, one would like to be able to compute 
the radiated noise field that takes into account 
these effects. One of the driving reasons be- 
hind this is the possibility of reducing the air- 
craft noise footprint by taking advantage of en- 
gine and wing location and manipulation of the 
nonuniform mean flow around the fuselage. In 
this paper, we present a spectral method for 
large scale computation of engine noise propa- 


gation in a realistic configuration and flow en- 
vironment around the aircraft. 

Recent advances both in computer archi- 
tecture and computing methods have already 
prompted researchers to address the fan tone 
noise problem using Computational Aeroacous- 
tics (CAA) time-domain techniques [10, 14, 15]. 
Stanescu et al. [15] used for this purpose a multi- 
domain spectral method and showed that ac- 
curate results can be obtained for noise prop- 
agation with only a relatively small (between 
four and five, depending on the highest polyno- 
mial degree) number of points per wavelength 
(PPW). The spectral multi-domain method is 
geometrically flexible, as it only requires finite- 
element type grids that can be generated for 
usual aircraft configurations with commercially 
available packages. Although this fact makes 
practical three-dimensional problems tractable, 
the computer resources needed for frequencies 
typical of those encountered in fan noise are still 
very large. Such problems can only be solved if 
a proper parallelization strategy is used in or- 
der to make best use of the modern distributed 
memory computers available. 

The Discontinuous Galerkin (DG) method 
for hyperbolic systems of conservation laws has 
been studied by Cockburn and Shu [2] in de- 
tail. A spectral element implementation was 
proposed by Kopriva et al. [8], and several stud- 
ies have been performed on its dissipation and 
dispersion properties [5, 11, 16]. These stud- 
ies indicate that it offers clear advantages in 
terms of accuracy when applied to wave prop- 
agation problems. The DG method is oth- 
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erwise completely similar to the multi-domain 
method [15] both as formulation and implemen- 
tation, and has therefore been considered as the 
discretization technique for this work. Through 
a proper choice of the flux computation points, 
the method only requires communication be- 
tween elements that have common faces, and 
is therefore very well suited for parallelization 
using message passing. In this setting, an ini- 
tial ordering of the elements pertaining to each 
processor allows communication of data for ele- 
ments on partition boundaries to be overlapped 
with computation of the fluxes for those ele- 
ments inside each partition, thus minimizing the 
impact of communication on the overall perfor- 
mance. 

The present paper is organized as follows. 
The next section briefly describes the DG 
method and its numerical implementation. Is- 
sues related to geometry representation, grid 
generation and parallelization strategy are all 
addressed. The numerical results section then 
presents computations of the sound field gen- 
erated by a single spinning mode propagat- 
ing around an axisymmetric nacelle, and com- 
pares the results with those obtained for prop- 
agation of the same mode when the nacelle is 
mounted in a generic fuselage-pylon-nacelle con- 
figuration, both with and without a wing. The 
conclusions section ends the paper. 


II. The DG Spectral Element Method 


1 Governing Equations 


The three-dimensional, nonlinear Euler equa- 
tions of gas dynamics in conservation form 


dQ A dFa 

dt h dx d 


= 0 


(1) 


are considered here to model the noise propaga- 
tion process. The state and flux vector compo- 
nents are, with the usual notation for density, 
internal energy, pressure and velocity compo- 


nents: 
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2 DG Discretization 


The computational domain is divided into non- 
overlapping general hexahedral elements that 
can have curved boundaries. Each element is 
then individually mapped onto the master el- 
ement II m = [ — 1 , l] 2 3 using an isoparametric 
transformation. Under the mapping, Eq. (1) 
becomes 


dQ dFd 


= 0 


(3) 


where the new variables are the transformed 
components of the state and flux vector 


Q = JQ , 

kJ •Xj'xxi 

m= 1 


(4) 


and J is the Jacobian of the transformation 


J = det 


9(xi,...,x 3 ) \ 
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(5) 


The computational space coordinates will be de- 
noted here by both (£i,£ 2 j£s) and for 

convenience. 


Let the space of polynomials of degree N in 
£ £ [—1,1] be denoted by Pm- A basis for this 
space can be constructed using the Lagrange in- 
terpolating polynomials 


hm 


n (*-fr) 
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( 6 ) 


through the N + 1 Gauss-Legendre [1] quadra- 
ture nodes i = 0,1,..., N. The three- 
dimensional solution Q is searched in a trial 
space obtained as the tensor product of one- 
dimensional degree N polynomials for each co- 
ordinate direction, P^ = Pn(0 x Pn(v ) x 
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(12) 


P/v(0) hence it takes the form 
N N N 

i=0j=0k=0 

(7) 

where Qijk(t) denote pointwise values of Q at 
time t. In view of a Galerkin discretization, the 
residual is required to be orthogonal to the same 
trial space within each element. Thus, the DG 
spectral element method looks for a solution of 
the system 


d^F= Y,Fi jk (Kihm)N 

m = 0 

where we introduced the discrete inner product 
N 

(u,v)n = y ^u n v n w n , the differential operator 

n = 0 

for £ for example can be written 

(13) 


{Qt, <t>ijk ) + (V^ • F, 4>ijk) — 0 (8) 

to be satisfied for all i,j, k = 0, 1, ... , N. Here 
(•, •) represents the usual L 2 inner product, and 
4>ijk = hi(^)hj(rj)hk(C) are the basis functions 
of Pjy. Using the divergence theorem, the equa- 
tion is recast as 


(.Qt:$ijk')F / (ftijkF'TldS — (P, ^ £$ijk) (9) 

J d^M 

In order to solve these equations efficiently, the 
integrals are computed using Gauss quadrature 



N 

f{O d t = 

m = 0 


V/ € P 2 N +1 


(10) 


w m being the Gauss weights. This replacement 
is exact for elements with straight edges, in 
which case J is a polynomial of degree one. In 
the general case of elements with curved bound- 
aries, a discretization error is incurred in the use 
of the numerical quadrature rule Eq. (10). 


Expanding the boundary integral and per- 
forming some algebraic manipulation, the final 
discrete form of the equations governing each 
variable at the Legendre- Gauss points can be 
found to be 


dQijk 

dt 


+ ZU + Dt 


F 


( 11 ) 


where the right-hand side is a sum of discrete 
differential operators acting on the flux values 
of an element, which include values on the ele- 
ment faces. Defining the contribution due to the 
volume term on the right-hand side of Eq. (9) 


and the remaining two follow by obvious per- 
mutations. In this equation the notation F* 
denotes a common face flux as explained below. 
The position of the Gauss-Legendre points in- 
side the element, as well as the position of the 
corresponding flux points on the element faces 
where F* must be evaluated, is shown in Fig. 1, 
drawn for one £ = const plane and N = 2. 

We note that the operation in Eq. (12) is 
actually a one-dinrensional operation that is ex- 
pressed in the computer code as the multiplica- 
tion of a matrix of entries = (L(, h m ) jy with 
a vector of size N + 1 containing the values of 
the F\ flux component at the Legendre-Gauss 
points on a line i = const inside the element. 
These flux values can be computed directly from 
the values of the state vector which are read- 
ily available. The values of the state vector at 
the face points of each element, i.e. Q{l,r)jXk) 
and Q(—l,r]jXk) for the £ = const faces, are 
obtained by evaluating the interpolating poly- 
nomial based on the Gauss-Legendre points at 
£ = ±1, again a one-dinrensional (implemented 
as a vector dot product) operation. However, 
these values can not be used directly in com- 
puting the pointwise flux values on the element 
faces, as they are not uniquely defined. Indeed, 
for each face point, two different values are ob- 
tained for the value of the state vector Q, one 
from each of the two neighboring elements. In 
order to define a unique ffux, a Rienrann prob- 
lem [13] is solved at each flux face point with 
the two different states as initial data, just as 
in the finite volume method. Letting the super- 
scripts L and R denote the two (left and right) 
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neighboring elements of a face, the application 
of the Riemann solver can be written 

FK-) = FKQ R (-),Q L (-)) (14) 

and completes the evaluation of the right-hand 
side of Eq. (11). Upon this, it becomes on ordi- 
nary differential equation which is integrated in 
time using a low-storage Runge-Kutta scheme 
optimized for wave propagation [17]. 

The computer code thus obtained can be 
used for computing both a steady mean flow, 
when appropriate inflow/outflow boundary con- 
ditions are specified at the fan face and the far 
field, and the propagation of acoustic waves su- 
perposed on such a mean flow previously ob- 
tained. For acoustic computations, once the 
mean flow is known, the boundary condition 
on acoustic source surfaces is allowed to be- 
come a function of time (corresponding to a 
source term in the governing equations), and 
the acoustic pressure is obtained as the differ- 
ence between the computed value of p and the 
mean flow pressure. Radiation boundaries are 
treated by a damping layer approach, similar 
to Stanescu [15]. The damping layer is acti- 
vated only for acoustic computations. This pre- 
cludes the need for a time-dependent radiation 
boundary condition, i.e. the boundary condi- 
tion at the exterior limit of the damping layer 
is the same inflow/outflow characteristic bound- 
ary condition that can be used to compute the 
steady mean flow field. 

3 Isoparametric Elements 

To obtain an accurate representation of the 
underlying geometry, the ICEMCFD Hexa [6] 
commercial package has been used for hexag- 
onal grid generation around the configurations 
of interest. Once an unstructured grid of hexa- 
hedra has been obtained, a new feature of this 
package allows for the generation of points along 
each of the edges of the hexahedral mesh ac- 
cording to a prescribed distribution, which can 
be either a Legendre-Lobatto or a Chebyshev- 
Lobatto distribution in arclength, for a specified 
polynomial degree N. The resulting points are 
written to a data file which can be read in a com- 


puter code and used to define the geometry with 
the same accuracy as the solution. All the nec- 
essary point coordinates can then be computed 
by interpolation based on the spectral inter- 
polants along the edges (to obtain coordinates 
at the Gauss points from the Lobatto points on 
the edges) followed by three-dimensional trans- 
finite interpolation [4] on the faces and inside 
the elements. The intermediary use of the Lo- 
batto points is found necessary because they in- 
clude the eight corners of the hexahedra (and 
hence the end points for each edge of the grid) 
which can thus be guaranteed to have a unique 
position. This would not be the case if their 
coordinates were created by interpolation from 
data at Gauss points (they would be different 
due to interpolation and truncation errors). 

4 Parallel Implementation 

The DG spectral element method is particu- 
larly well suited for a parallel implementation 
through message passing as the computation is 
naturally structured by elements, with most of 
the work done inside the element i.e. evalua- 
tion of derivatives in Eq. (11). Here we briefly 
describe our parallel implementation. The pro- 
gram starts with a preprocessing step wherein 
the underlying hexahedral mesh is decomposed 
in as many partitions as the number of pro- 
cessors using the multilevel graph partition- 
ing techniques in METIS [7]. Each processor 
creates the connectivity for its grid partition, 
which includes a list of elements (’’cut elements” 
henceforth) having one or more ’’cut faces” in 
common with elements in neighboring parti- 
tions. 

One stage of the Runge-Kutta time integra- 
tion scheme for the parallel DG algorithm can 
then be described as follows: 

• Compute state vector values at the face 
flux points through interpolation for all 
cut elements. 

• Post non-blocking send/receive for the 
state vector values on the cut faces. Each 
processor sends the values interpolated 
from the cut elements it owns to the cor- 
responding neighboring partition. 
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• Compute state vector values at the face 
flux points through interpolation in the 
regular elements inside each partition. 

• Compute unique fluxes through the Rie- 
rnann solver for regular faces, Eq. (14). 

• Test for completed send/receive. 

• Compute unique face fluxes through the 
Riemann solver for cut faces, Eq. (14). 

• Compute residual and update the state 
vector according to Eq. (11). 

The parallel performance of this algorithm 
has been tested on several architectures, for var- 
ious grid sizes and polynomial degrees N. We 
present here the results obtained for two grid 
sizes. The first one is a small grid of 884 ele- 
ments with polynomial degree N = 12, which 
has been run on an IBM SP3 machine with 
375-MHz Power3 processors, each with 512MB 
of memory. The second one is a medium-size 
(17370 elements) grid used for the fuselage- 
nacelle configuration, run with N = 4 on an 
SGI 0rigin2000 machine with 250MHz R10000 
processors. Both cases could be run on a sin- 
gle processor, and Fig. 2 shows the performance 
on up to 16 processors as obtained by compar- 
ing the wall-clock time with the single processor 
run. 

III. Numerical Examples 

The methodology described above has been val- 
idated by comparison with experiment in [15]. 
The computations presented here have been 
performed in order to study the parallel per- 
formance of the method, the problems raised by 
such large scale computations, as well as to gain 
some initial physical insight into the influence of 
the fuselage and/or other surfaces on the sound 
field of a nacelle. Only inlet noise is modeled 
in all cases. To keep the paper at a reasonable 
length, results are presented only for a station- 
ary medium. Mean flow effects will be discussed 
in another work. 

Remarks on non-dimensionalization. In all 
cases, lengths were non-dimensionalized using 


the radius of the inlet duct as reference value. 
Sound pressure levels have been computed as 
SPL = 201og 10 (p RM s) + 100 where jjrms is the 
root-mean square pressure non-dimensionalized 
by PooC^o- This is appropriate since the com- 
putations for these configurations were mainly 
intended to prove the feasibility of the approach, 
not to compute the absolute levels of noise. The 
same amplitude of the incoming acoustic wave 
was used for all computations. In the figures, 
SPL contours in the range —20 dB through 
20 dB are plotted with an increment of 5 dB. 

1 Nacelle Alone 

The first computation was performed for the 
sound field radiated from an axisymmetric bell- 
mouth nacelle. The radius of the leading edge 
of the bell-mouth has been chosen to be 1/4 of 
the radius of the inlet duct, and no center-body 
(hub) was added. Radiation of the (2, 1) spin- 
ning mode at reduced frequency u> = 8.0 in a 
quiescent fluid has been modeled on a grid made 
of E = 9919 elements. The polynomial degree 
used was N = 6, leading to an average PPW 
value of about 5 (largest element edge length 
is 1.06). There are a total of 3,402,217 Gauss- 
Legendre points in the domain. The nacelle cas- 
ing extends on the x-axis (the symmetry axis) 
from x = —1.25 to x = 3.0, and the domain 
extends from x = — 15 to x = 8, with a damp- 
ing layer width of 3 in each direction. Incoming 
acoustic waves are specified on the source sur- 
face which is the circular disk with center at 
(0,0,0) and radius 1. 

Fig. 3 shows the SPL contours on the surface 
of the nacelle and in a plane that cuts the na- 
celle through its symmetry axis. Since no scat- 
tering surfaces other than the engine exist, the 
radiation is symmetric with respect to the en- 
gine axis. Fig. 4 shows the contours in the plane 
z = —4. This computation was also used as a 
test of the radiation boundary conditions, which 
perform quite well in this case (very little wavi- 
ness of the SPL contours is noticeable towards 
the interior limits of the damping layer). The 
main radiation lobe makes an angle of about 
9 = 46 degrees with the x-axis, a value signif- 
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icantly different from 8 = 57 degrees obtained 
using the analytical formula derived by Rice [12] 
for cylindrical ducts, 

0088 = ^1-^- (15) 

where k mfJj = 6.70613 is the circular duct eigen- 
value for this mode. This difference is most 
probably caused by the large thickness of the 
nacelle lip (half of the inlet duct radius). 

2 Fuselage-Nacelle Configuration 

The second case considered a configuration ob- 
tained by joining the same nacelle with a cylin- 
drical fuselage through an elliptic cylinder py- 
lon. The cylinder representing the fuselage has 
a radius of 4.5 and its center is situated at 
y = — 7. The grid is made up of E = 17370 
elements of the same order N = 6 (5,957,910 
Gauss-Legendre points). The minimum and 
maximum values of the coordinates are (-18,- 
7,-9) and (8,10,8) respectively, and the damp- 
ing layer has the same width as in the previous 
case, with the difference that the plane y = —7 
does not have an adjacent damping layer in the 
direction perpendicular on it. Instead, a zero 
normal acoustic velocity boundary condition is 
used on this plane to simulate the effect of a 
symmetric configuration with engines on both 
parts of the fuselage. 

As can be seen in Fig. 5, the presence of the 
fuselage creates a series of secondary radiation 
lobes due to the reflection of the nacelle lobes 
from the surface and to the various interferences 
that occur. The structure of the radiated sound 
field becomes thus quite complex. The presence 
and extent of the footprint of the SPL = 0 con- 
tour level in this plane is a measure of the effect 
of the fuselage in increasing the noise levels com- 
pared to the nacelle alone case. The figure also 
shows the trace of mesh decomposition on the 
solid surfaces defining the geometry. 

3 Fuselage- Wing-Nacelle Configuration 

The geometry used for the previous case was 
augmented with a wing of elliptic cross-section, 


slightly swept and mounted below the nacelle. 
The center of the ellipse defining the cross- 
section of the wing is situated in the plane 
z = —3. A grid of E = 22843 elements was gen- 
erated in this case (for N = 6 this corresponds 
to 7,835,149 discretization points). In all other 
respects, the definition of the domain is simi- 
lar to the previous case. Shown in Fig. 6 are 
the SPL contours on the solid surfaces defining 
the geometry for this case. The presence of the 
wing is responsible for a complex, considerably 
asymmetric, sound field due to sound scattering 
from its leading and trailing edges. As can be 
expected however, although in the plane z = 0 
above the wing the noise levels are higher than 
in the preceding case, in a plane below the wing 
the noise levels are remarkably lower, due to the 
upward reflection of the sound field by the wing 
surface, the difference immediately below the 
wing lying in the 10 — 15dB range. The corre- 
sponding plots are shown in Figs 7 and 8. Also 
obvious from these figures is that the damping 
layer region must be extended for a more ac- 
curate computation, as the contours are jagged 
towards the limits of the domain. 


IV. Conclusions 


The qualitative results presented here show that 
a more complete and accurate modeling of tone 
noise radiated from aircraft engines, including 
the effect of scattering from wing and fuselage, 
is now possible by use of the large distributed 
memory parallel computers available nowadays. 
This can have major implications in the de- 
sign stage, in terms of affecting the noise foot- 
print of the aircraft by changing the relative 
engine/wing positioning and altering the wing 
planform. The paper discusses a time domain 
methodology that is being developed to perform 
such studies. Further work, including the effect 
of mean flow and the sound radiated from the 
bypass exhaust, will be presented in the near 
future. 
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Fig. 1: A £ = const plane through the master element, showing the position of the Gauss-Legendre points 
inside the element and the position of the flux points on the element faces. 



Fig. 2: Parallel performance for two test cases on up to 16 processors. 
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Fig. 3: SPL contours for radiation of spinning mode (2, 1) at u> = 8 from an axisymmetric bell-mouth nacelle 
alone. Nacelle surface and plane z = 0 through its axis shown. 



Fig. 4: SPL contours for radiation of spinning mode (2, 1) at to = 8 from an axisymmetric bell-mouth nacelle 
alone. Plane 2 = —4 shown. Trace of mesh decomposition shown on nacelle surface. 
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Fig. 5: SPL contours for the radiation of spinning mode (2,1) at u = 8 in the fuselage-pylon-nacelle 
configuration in plane z = — 4 . Trace of mesh decomposition shown on the solid surfaces. 



Fig. 6: SPL contours for the radiation of spinning mode (2,1) at u> = 8 in the complete configuration with 
wing. Solid surfaces shown only. 
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Fig. 7: SPL contours for the radiation of spinning mode (2,1) at u> = 8 in the complete configuration with 
wing. Plane z = 0 and solid surfaces shown. 



Fig. 8: SPL contours for the radiation of spinning mode (2,1) at ui = 8 in the complete configuration with 
wing. Plane 2 = —4 and trace of mesh decomposition on solid surfaces shown. 
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